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Abstract: We extract at redshift z = a. Milky Way sized object including gas, stars 
and dark matter (DM) from a recent, high-resolution cosmological N-body simulation with 
baryons. Its resolution is sufficient to witness the formation of a rotating disk and bulge 
at the center of the halo potential, therefore providing a realistic description of the birth 
and the evolution of galactic structures in the ACDM cosmology paradigm. 
The phase-space structure of the central galaxy reveals that, throughout a thick region, 
the dark halo is co-rotating on average with the stellar disk. At the Earth's location, the 
rotating component, sometimes called dark disk in the literature, is characterized by a min- 
imum lag velocity viag — 75 km/s, in which case it contributes to around 25% of the total 
DM local density, whose value is pdm — 0.37 GeV/cm^. The velocity distributions also 
show strong deviations from pure Gaussian and Maxwellian distributions, with a sharper 
drop of the high velocity tail. 

We give a detailed study of the impact of these features on the predictions for DM signals 
in direct detection experiments. In particular, the question of whether the modulation 
signal observed by DAMA is or is not excluded by limits set by other experiments (CDMS, 
XENON and CRESST...) is re-analyzed and compared to the case of a standard Maxwellian 
halo. We consider spin-independent interactions for both the elastic and the inelastic 
scattering scenarios. For the first time, we calculate the allowed regions for DAMA and 
the exclusion limits of other null experiments directly from the velocity distributions found 
in the simulation. We then compare these results with the predictions of various analytical 
distributions. 

We find that the compatibility between DAMA and the other experiments is improved. 
In the elastic scenario, the DAMA modulation signal is slightly enhanced in the so-called 
channeling region, as a result of several effects that include a departure from a Maxwellian 
distribution and anisotropies in the velocity dispersions due to the dark disk. For the 
inelastic scenario, the improvement of the fit is mainly attributable to the departure 
from a Maxwellian distribution at high velocity. It is correctly modeled by a general- 
ized Maxwellian distribution with a parameter a ~ 1.95, or by a Tsallis distribution with 
q ~ 0.75. 

Keywords: dark matter, N-body simulations, direct detection experiments. 
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1. Introduction 

Nowadays, the hypothesis of dark matter (DM) is considered as the most plausible expla- 
nation for various observations from galactic to cosmological scales. Flat rotation curves in 
spiral galaxies [1,2], velocity dispersions in galaxy clusters [3,4], gravitational lensing mass 
reconstructions [5-7], hierarchical patterns in large scale structures [8,9] and cosmic mi- 
crowave background anisotropics [10] all point towards the existence of a new neutral, non 
baryonic, cold (i.e. with low thermal velocities) and weakly interacting massive particle 
(WIMP) in the universe. 

Theories beyond the Standard Model of particle physics [11-13] provide many candi- 
dates for DM particles (see for example Refs. [14-17]). However, its existence will not be 
proved and its properties will not be determined until a clear identification is made through 
signals in indirect or direct detection experiments or in particle accelerators. Recently, sev- 
eral astrophysical excesses have been reported [18-21], and associated with a DM signal 
in more or less exotic scenarios (see, among many others, Refs. [22-27]). However, up to 
now, none of them requires DM as an irrefutable explanation. The INTEGRAL 511 keV 
line 7-ray signal [18] can be explained by the emission of low mass X-ray binaries [28]. 
The EGRET diffuse 7-ray GeV anomaly [19] was most likely due to an error in the en- 
ergy calibration [29], and has since indeed been ruled out by the FERMI instrument [30]. 
The positron and electron peak seen by ATIC [21] has been flattened by HESS [31] and 
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FERMI [32] measurements. The positron fraction rise observed by PAMELA [20], although 
striking, can be accommodated in standard astrophysical scenarios with pulsars [33,34], or 
with an inhomogeneous distribution of supernovae remnants [35]. 

Direct detection provides a unique opportunity to disentangle DM signals from other 
astrophysical processes. The idea is to directly measure the energy deposited in a low 
noise detector during the collision between a DM particle and a detector nucleus. Several 
experiments have been designed for this aim [36-40]. Up to now, all of them are compatible 
with null results, except the DAMA/LIBRA collaboration which claims to have observed 
a signal with an annual modulation characteristic of an DM signal [41]. This modulation 
is induced by the Earth rotation around the Sun, and is presented as a signature hard to 
reproduce by background events. The DAMA result is still quite controversial because it 
leads to tensions with other experiment measurements (for a recent status, see Refs. [42- 
45]). 

Predicting direct detection signals in any particle physics model requires astrophysical 
assumptions in the form of a value for the local DM density as well as velocity distributions 
in the solar neighborhood. Usually, the local DM density is taken as poM = 0.3 GeV/cm^ = 
0.079 Msun/pc'^, a value which is supported by the Milky- Way rotation curve [46,47], and 
stars spectrophotomctric data [48]. A recent determination which takes numerous dynam- 
ical obscrvablcs into account gives a larger value, poM = 0.39 GcV/cm'^ [49]. However, it 
should be mentioned that determinations based on mass models only give average values. 
A local peak in the DM density, although improbable without a dynamical reason [50], 
is in principle not excluded by very local bounds based on planetary data [51,52]. An- 
other common over-simplified assumption is to take the velocity distribution as isotropic 
and Maxwellian. Incomplete virialization, streams and interaction with the stellar disk 
could contribute to create anisotropies or a non Maxwellian spectrum. Moreover, even 
in a completely virialized structure at equilibrium, anisotropy is expected as a result of 
the inhomogeneous gravitational potential with a power law density profile [53,54]. The 
consequences for direct detection of some of these departures from the standard halo case 
have been analyzed by several authors [55-62]. 

Numerical simulations have been used for decades as a powerful virtual laboratory to 
explore complex dynamical systems. On cosmological scales, early studies enabled to assess 
the role of DM in the hierarchical structure formation sheme [63,64]. Recent techniques [65, 
66] as well as improving computing capabilities have brought the possibility to extract 
galactic DM halos with tremendous resolution. The so-called Via Lactea simulation [67] 
or the Aquarius project [68] have resolved a Milky- Way sized galactic halo with more 
than a billion particles. It appears that numerous sub-halos are present, anisotropies and 
deviations from standard Maxwellian velocity distributions are also found and lead to 
somewhat different signatures in direct detection experiments [43,53,54,62,69]. Despite 
their impressive resolution, these DM only simulations fail as a realistic description of a 
galactic halo for the simple reason that they totally neglect the baryonic components (stars 
and gas in the galactic disk and bulge), although these dominate are known to dominate 
the dynamics near the galactic center. This shortcoming is being fixed by very recent 
cosmological simulations which include baryons [70,71]. It appears that the DM spatial 
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distribution is neither spherical, or triaxial, but has a thick oblate shape around the galactic 
disk, in what some authors have called a dark disk [72,73], that is co-rotating with the 
galactic stellar disk. These particular characteristics lead to a potential enhancement of 
the direct detection signal at low recoil energies [73] , although there is an important spread 
in the predictions, depending on the particular merger history for each galaxy. 

In this paper, a recent and advanced cosmological simulation with baryons is used as a 
realistic framework to extract detailed predictions for direct detection. The paper is orga- 
nized as follows. In Sec. 2, the simulation is presented and described. Velocity distributions 
in the solar neighborhood are extracted, analyzed and compared with standard simplified 
assumptions. The presence of a co-rotating dark disk is also discussed. In Sec. 3, after a 
brief overview of the event rate formalism and the experimental status, the direct detection 
predictions are given for both the elastic and the inelastic scenarios. Finally, our results 
are summarized in Sec. 4. 

2. A cosmological simulation with baryons 

Although galaxy formation is far from being completely understood, both theory and ob- 
servation have made significant progresses in the recent years. Based on first principles, it 
is now possible to perform simulations of Milky- Way-like galaxies in reasonable agreement 
with observed galaxies of similar circular velocities [74-76]. Some long-lasting problems 
seem however to remain, especially for massive galaxies like our own: the simulated bulge 
seems systematically too massive (sometimes referred to as the angular momentum prob- 
lem) , and the total baryon fraction too large (the overcooling problem) . Following previous 
work initiated by Bruch et al. [73], we nevertheless would like to use a Milky- Way-like 
galaxy simulation including both dark matter and baryons dynamics to infer the expected 
signal in dark matter direct detection experiments. 

2.1 Features 

The simulation was performed using the cosmological Adaptive Mesh Refinement (AMR) 
code RAMSES [77]. We follow the evolution of 2 collisionless fiuids, namely dark matter 
and stars, coupled through gravity to the dissipative baryonic component. Gas dynamics is 
described using a second-order unsplit Godunov scheme. Gas is allowed to cool radiatively, 
based on a standard equilibrium photo-ionised mixture of Hydrogen and Helium [78], taken 
into account the excess cooling due to metals [79]. Star formation is implemented using a 
standard Schmidt law, spawning stars as a Poisson random process [80]. Our rather stan- 
dard recipe for supernovae feedback and the associated chemical enrichment is described 
in Ref. [81]. 

We used the so-called "zoom-in" simulation technique, for which we identify a candi- 
date halo using a low resolution, dark matter only simulation, and then re-simulate it at 
higher resolution, degrading both spatial and mass resolutions as we move further away 
from the central region. We checked that our final halo is not contaminated by high mass, 
low resolution dark matter particles. Our eff'ective resolution in the initial conditions cor- 
responds to a Cartesian grid of 1024^ elements covering our 20 Mpc//i periodic box. The 
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Figure 1: View of the simulation box (top), with a zoom on the central and better resolved region 
(bottom). The size of the simulation box is 20 Mpc/h. The central galaxy is embedded in a large 
scale filamentary structure. 



corresponding dark matter particle mass is therefore 7.46 x 10^ solar masses, given we used 
the following cosmological model = 0.3, = 0.7, fi^ = 0.045, ffo = 70 km/s/Mpc. 
The initial Gaussian random field was generated according to the Eisenstein and Hu (1999) 
ACDM model [82] for the transfer function. 

At redshift z = 0, the complete simulation volume contains 8,376,152 dark matter 
particles of various masses, 11,702,370 star particles and 19,677,225 AMR cells. We have 
fixed the spatial resolution to be 200 pc in physical units, which ensures that the vertical 
scale height of the stellar disc is barely resolved. The final halo Virial radius Rvir^ is 
found to be 264 kpc. It contains Nqm = 842, 768 dark matter particles corresponding 

^The Virial radius Rvir is defined as a function of the cosmological matter density pm by '^'^^^'^^i"^ = 
2mpm 
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Figure 2: View of the dark matter halo of the central galaxy of Fig. 1. This halo will he used in 
all analyses in this paper. The radius of the red circle is equal to its virial radius, i.e. 264 kpc. As 
expected, this halo is centrally concentrated and has a large number of subhalos. 



to a halo mass Mh = 6.3 lO^Msun- Within the central 60 kpc, we have A'^ = 5, 416, 865 
stars particles (within 15 kpc, this number goes down to A^* =4,927,325 particles). We have 
carefully selected our halo so that its mass accretion history is typical of an late-type galaxy, 
with no major merger and a steady accretion rate in the last 8 Gyr. Moreover, our halo 
mass is small enough (6 x 10^^ M0) so that its initial Lagrangian volume, with a comoving 
radius of 2 Mpc/h, is not affected by spurious effects due to the periodic replication of the 
20 Mpc/h box. Our simulated halo, as can be seen in Figure 1, is embedded in a large 
scale filament, which is quite typical for a galaxy. 

As discussed above, our central galaxy is slightly too concentrated, with a bulge mass 
of the order of 4xlO^'^Msun and a similar disc mass, while the most recent observational 
constraints in the Milky Way suggest a bulge mass of 2xlO"'^'^Msun for a disc mass of 
6xl0i°Msun [47]. Although the total mass in stars of our simulated galaxy agrees quite 
well with observations, we stress that we have considered here a rather low mass halo 
model {Myir ~ 6 x lO^^Msun), resulting in a total baryon fraction in the galaxy close to 
the universal value of 15%. One should bear in mind that the larger bulge mass of our 
simulated galaxy will have the effect of concentrating the dark matter distribution further 
in. This process, usually referred to as "Adiabatic Contraction" (AC), results in a decrease 
of the dark matter characteristic scale radius, and in an increase of the dark matter circular 
velocity [83-85]. When the effects of AC are combined with a small total halo mass (6 x 10^^ 
Mq), we obtain a dark matter density in the solar neighbourhood close to 0.37 GeV/cm'^. 
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c) Gas face on d) Gas edge on 



Figure 3: Face-on and edge-on distributions of the baryonic material in our galaxy. The white 
vertical line in the upper left panel has a length of 10 kpc. Note that both the stars (upper panels) 
and the gas (lower panels) have formed a thin disk, in good agreement with observations. A number 
of satellites can also be seen, more numerous in the stars than in the gas, the gas having been 
expelled from a large fraction of the substructures. The filamentary structures in the gas edge- on 
view outlines past satellite trajectories. 



Reducing the bulge mass by a factor of 2 would result in a 10% decrease of the circular 
velocity and in a 20% decrease of the dark matter density in the solar neighbourhood (see 
Ref. [84] for a quantitative estimate). 

The central dark halo, stars and gas are shown respectively on figure 2 and 3. Figure 2 
shows the dark matter distribution. As expected, it is centrally concentrated and has a 
considerable number of subhaloes. Using Adaptahop to extract the substructure informa- 
tion, we find that the central galaxy contains 92 subhaloes, of which 79 are inside the virial 
radius. The subhalo mass fraction is 4.8%, a rather small value due to the limited mass res- 
olution in the simulation. The closest clump to the galactic center is at a radial distance of 
12.7 kpc. Therefore, as no clump lies in the the solar neighborhood region (7 < -R < 9 kpc 
and — 1 < z < 1 kpc), the influence of substructures on DM direct detection signals cannot 
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be inferred from this simulation. 

Figure 3 shows the distribution of the stars and the gas. Both display a very thin disk 
of a radial extent of roughly 10 kpc and a clear spiral structure. For the stellar population, 
this is mainly seen in the outer parts of the disc, while for the gas it can be traced to 
a region very near the center. In both cases it is two-armed, but has clear asymmetries. 
Secondary branchings of the arms, coupled to the main arms result in a ring-like shape in 
the gas distribution. 

Note also that there is an important central concentration in both distributions. This 
is of short extent, both in the equatorial plane and perpendicular to it. The latter property, 
as well its rich gaseous content which has roughly the same extent as the stellar distribution, 
exclude the possibility that it is a classical bulge and argue strongly that is a discy bulge [86] . 
Further investigations of this point is nevertheless necessary. 

In the stellar distribution the main galaxy is surrounded by a number of small conden- 
sations, big enough to be considered as small satellite galaxies. In the gaseous distribution, 
however, only one of these condensations can be seen. This means that the remaining 
condensations have no gaseous counterpart; the gas having been expelled either by ram 
pressure stripping, or by a central supernova, or a combination of both. The clear filamen- 
tary structures in the gaseous edge-on plot outline past satellite trajectories. 

2.2 Velocity distributions 

In this section, we extract and analyze the velocity distributions given by our N-body 
simulation, and compare them with the usual Gaussian/Maxwellian assumption. As the 
goal of this study is to calculate direct detection predictions with a realistic DM halo, the 
emphasis will be put on the local phase-space structure, in the solar neighborhood. 

The Sun being located at a distance i?o — 8 kpc from the galactic center, we first select 
all DM particles in the spherical shell 7 kpc < R < 9 kpc. The number of particles in the 
selection is Ngheii = 16,545. Following Ref. [69], we calculate the velocity distributions 
along the principal axes of the velocity dispersion tensor afj. They are shown on Fig. 4. 
The Gaussian fit for each component and the Maxwellian fit for the velocity module are 
shown in red. A clear departure from these distributions is seen. 

In order to model and parameterize these deviations, one can introduce generalized 
Gaussian and Maxwellian distributions, which are convenient fitting functions. The gen- 
eralized Gaussian distribution with a mean velocity //, a velocity dispersion parameter vq, 
and a Gaussianity parameter a is defined by 



with a normalization factor N(vo,a) = 2voT{l + l/2a). For a = 1, the usual Gaussian 
distribution is recovered. Similarly, the generalized Maxwellian distribution can be defined 



with a normalization factor N{vo,a) = 47rT;Qr(l + 3/2a)/3. For a = 1, we get the usual 
Maxwellian distribution. 





by 
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Figure 4: Velocity distributions of dark matter particles (Nsheii = 16, 545 ) in a spherical shell 
7 < R < 9 kpc around the galactic center. 

a), b) and c) Velocity components along the principal axes of the velocity dispersion tensor, together 
with the Gaussian (red) and a generalized Gaussian (green) distribution fits (cfr. Eq. (2.1)). 
d) Velocity module, with Maxwellian (red), Tsallis (green) and generalized Maxwellian (orange and 
purple) fits (cfr. Eqs. (2.2,2.3)). 

jjt, a (both in km/s) and K stand for the mean, the standard deviation and the Kurtosis parameter 
of the distribution. The goodness of fit is indicated by the value of the ^ vs. the number of degrees 
of freedom ( dof ) . 



For the velocity components, a fit with a generahzed Gaussian distribution (in green 
in Fig. 4, panels a), b) and c)) shows that the velocity distribution given by the simulation 
is systematically more flat than a Gaussian distribution. This property can be described 
by the so-called Kurtosis parameter which compares the fourth-order moment with the 
square of the variance. A Kurtosis parameter < 3 corresponds to a distribution that is 
platykurtic, i.e. more flat than a Gaussian distribution with the same standard deviation. 
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Figure 5: Velocity distributions of dark matter particles (Nring — 2, 662J in a ring 7 < R < 9 kpc, 
\z\ < 1 kpc around the galactic plane. 

a) Radial velocity Vr, with Gaussian (red) and generalized Gaussian (green) fits (cfr. Eq. (2.1)). 

b) Tangential velocity v^, with a double Gaussian fit. f indicates the fraction of each component. 

c) Velocity across the galactic plane Vz, with Gaussian (red) and generalized Gaussian (green) fits 
(cfr. Eq. (2.1)). 

d) Velocity module, with Maxwellian (red) and a generalized Maxwellian (green) fit (cfr. Eq. (2.2)). 
fi, a (both in km/s) and K stand for the mean, the standard deviation and the Kurtosis parameter 
of the distribution. The goodness of fit is indicated by the value of the 'vs. the number of degrees 
of freedom ( dof). 



while K > 3 corresponds to a distribution that is leptokurtic, i.e. more peaked around the 
central value compared to the Gaussian one. So velocity distributions extracted from this 
simulation are systematically platykurtic, in contrast with results obtained in DM only 
simulations [69,87]. 

The goodness-of-fit for each theoretical distribution is calculated using the Pearson's 
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test. The values of (calculated with a number of bins Uhin = 50) printed in each 
panel of Fig. 4 show that the fit is considerably improved when using a generalized Gaus- 
sian distribution with a > 1 (the actual values of vq and a are estimated on the basis of 
the standard deviation a and the Kurtosis parameter K extracted from the velocity distri- 
bution) . We also notice that along the principal direction "3" , the velocity dispersion and 
the deviation from the Gaussian distribution are smallest. It appears that this direction is 
actually very close to the galactic pole axis. 

For the velocity module, a clear deviation from a Maxwellian distribution (in red) is 
apparent, with a sharper drop of the high velocity tail, particularly crucial for the inelastic 
recoil scenario (see Sec. 3.5). There is no evidence of bumps in the velocity distribution, as 
seen in Ref . [69] . A fit of the distribution with a generalized Maxwellian distribution with 
a > 1 leads to some improvement of the goodness-of-fit. However, no value of a enables 
to model the simulation distribution in a satisfactory way. With two free parameters, 
only the mean value and the standard deviation can be accommodated, while the Kurtosis 
parameter is derived. It appears that the N-body velocity distribution is "fatter" than the 
best-fit generalized Maxwellian (vq = 301.6 km/s and a = 1.55), with an observed Kurtosis 
parameter K = 2.39 compared to K = 2.71 for the fit. In panel d) of Fig. 4, we show 
two possible fits with a generalized Maxwellian distribution, corresponding to a = 1.5 and 
a = 1.95. The former gives a better global fit for the whole distribution and for the low 
velocity bins, the latter gives a better fit of the high velocity tail. 

In fact, the sharp drop of the high velocity tail has been recognized as a universal 
behavior of relaxed collisionless structures [54]. The presence of long-range gravitational 
forces in dark matter structures indicates that these systems should be described by non- 
extensive statistical mechanics. The generalization of the usual Boltzmann-Gibbs entropy 
to non-extensive systems by Tsallis [88] leads to distribution functions of the form 



where N(vQ,q) is a normalization constant. The Maxwell-Botzmann distribution is recov- 
ered by taking the limit g — )• 1. Equilibrated self-gravitating collisionless structures have 
been shown to exhibit Tsallis distributions both analytically [89], and numerically [54,90]. 
For the particles in the spherical shell around 8 kpc in this simulation, the velocity module 
distribution is indeed very well fit by a Tsallis distribution with vq = 267.2 km/s and 
q = 0.773. With these values, the Kurtosis parameter for the distribution is K = 2.44, 
very close to the observed value. 

As argued in Ref. [53], there should be some anisotropy in the velocity distributions 
between the radial and the tangential directions, as particles moving in the radial direction 
see a change in the potential. The velocity (dispersion) anisotropy is usually described 
by the parameter /3 = 1 — crf/cr^. Also, the departure from a Maxwellian distribution, 
parameterized by 1 — g for Tsallis distributions could show the same anisotropy. In this 
simulation, for the DM particles in the spherical shell 7 < i? < 9 kpc, we find a velocity 
anisotropy (3 ~ 0.12, together with a slightly smaller Kurtosis parameter in the radial 
direction K = 2.41, in general agreement with the conclusions drawn from the universal 
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Figure 6: Velocity distributions of star particles (Ngtar — 143, 320J in a ring 7 < R < 9 kpc, 
\z\ < 3 kpc around the galactic plane. 

a) Radial velocity Vr with Gaussian (red) and generalized Gaussian (green) fits (cfr. Eq. (2.1)). 
h) Tangential velocity Vrf, with a double Gaussian fit, with a thin (red) and a thick (green) disk 
components, f indicates the fraction for each component. 

fx, a (both in km/s) and K stand for the mean, the standard deviation and the Kurtosis parameter 
of the distribution. The goodness of fit is indicated by the value of the I'S- the number of degrees 
of freedom (dof). 



relation between the density slope and the velocity anisotropy presented in Ref. [53]. To 
calculate the velocity anisotropy, the radial velocity dispersion ar is directly extracted from 
the radial velocity distribution, while the tangential velocity dispersion at can be deduced 
from the mean ^ and the standard deviation a of the distribution of the velocity module 
as fj^ + 2(7^ = /i^ + cr^. 

In dark matter-only simulations, one has a 47r uncertainty on the Sun position. This 
freedom is lifted, at least partially, in simulations with baryons where the galactic plane 
can be identified. To determine it, we select star particles with 3 < ii < 10 kpc to avoid 
the galactic bulge. We have Ngtar = 1) 935, 104 in this selection. The galactic plane is then 
rotated onto the xy plane by diagonalizing the position tensor. In this new reference frame, 
denoted by {1^, ly, 1^}, we select DM particles with 7 < i? < 9 kpc and \z\ < 1 kpc, corre- 
sponding to particles located in a ring around the galactic plane. The number of particles 
in this selection is Nring = 2,662. Fig. 5 shows the corresponding velocity distributions 
in cylindrical coordinates {r, z}. Anisotropy is manifest. There is a significant rotation 
along the galactic disk, with (f,^) = 35.2 km/s. For the radial and z velocity components, 
the mean of the distribution is compatible with zero. The velocity dispersions and Kurtosis 
parameters also exhibit some anisotropy in the z direction compared to the galactic plane. 
For the distribution of the velocity module, a strong deviation from a Maxwellian distri- 
bution is again visible, with a cut-off of the high velocity tail. However, despite a coarser 
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Figure 7: (Left) Average density of DM particles with 7 < R < 9 kpc as a function of the height 
from the galactic disk z (R is the spherical radius to the galactic center). The dashed line gives 
the average value for the entire spherical shell. To select particles in z slices, we used a thickness 
6z — 2 kpc. 

(Right) Ratio of ring to shell densities as a function of distance from the galactic center for different 
planes. The ratio fluctuates around 1.2 for the galactic plane (blue), while it drops to a value ~ 0.9 
for other planes (green, magenta). For the plane yz, the sudden peak ai i? ~ 13 kpc is due to the 
presence of a satellite halo, visible on Fig. 8.b. 

resolution, we can observe that the distribution for ring particles differs significantly from 
the Tsallis form that is expected for equilibrated systems. Therefore, we can already stress 
that predictions of direct detection signals that use equilibrium distributions should be 
taken with some caution. 

2.3 Dark disk, halo rotation and local density 

As pointed out in Ref. [72], the rotation of the DM halo observed in Fig. 5 is expected in 
ACDM cosmology with hierarchical structure formation. At relatively low redshift, z < 2, 
merging satellite galaxies are preferentially dragged towards the already well-formed disk 
of a host galaxy, where they are disrupted by tidal forces. The accreted material settles 
in a thick stellar disk and a thick dark matter disk, dubbed the dark disk. The vertical 
distribution of the stars of the host galaxy thickens due to the impact, but all new stars 
form in a thin disc [91-93]. Observational evidence for a thick stellar disk in the Milky 
Way has been first pointed out by Gilmore and Reid (1983) [94]. Although the detailed 
morphology and some quantitative properties of the final galaxy depend on its merger 
history, the existence of a thick stellar disk and a co-rotating dark disk in the ACDM 
paradigm seems quite robust. Using three cosmological hydrodynamics simulations. Read 
et al. [71] conclude that the accreted dark matter can contribute ~ 0.25 — 1.5 times the 
non-rotating halo density at the solar position. Moreover, the DM distribution is best-fit 
with a double Gaussian, one for the non-rotating halo component, and one for the rotating 
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Figure 8: Density maps of the dark matter halo in the planes a) xy (galactic plane), b) yz. 
Contours correspond to pdm = {0.1,0.3, 1.0,3.0} GeV/cm'^. 
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Figure 9: Mean tangential velocity (v^) as a function of R for 2 = (left panel), and as a function 
of z for r ~ 8 kpc (right), where r — y'' + is the polar radius. 



accreted DM. The rotation lag of the rotating component is in the range — 150 km/s. A 
large value of the rotation lag, corresponding to a small halo circular speed, is only found 
for galaxies which had no significant merger after a redshift z = 2. The importance of 
the dark disk for the prediction of DM direct detection rates has been acknowledged by 
several authors [71,73]. It could lead to an increase by up to a factor of 3 in the 5 — 20 keV 
recoil energy range. The signal modulation can also be boosted and the modulation phase 
is shifted. However, in Ref. [95], the authors use high-resolution simulations of accretion 
events in order to bracket the range of co- rotating accreted dark matter. They find that 
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the Milky Way's merger history must have been unusually quiescent compared to median 
ACDM expectations. As a result, the fraction of accreted dark matter near the Sun is 
less than 20% of the host halo density, and does not change the likelihood of detection 
significantly. 

To investigate this topic further, we select star particles with 7 < i? < 9 kpc and 
\z\ < 3 kpc, we have Ngtar = 143,320. The distribution of Vr and are shown on 
Fig. 6. We observe that the dark matter and the star particles are indeed co-rotating 
in the solar neighborhood. The mean tangential velocity is (v^) = 201 km/s but tends 
towards (u^) = 225 km/s for stars closer to the galactic plane, which is consistent with 
Milky Way rotation curve data [47]. Moreover, the distribution is clearly bimodal, a 
double Gaussian gives a reasonable fit. The low velocity component (/i = 175 km/s) has 
a larger velocity dispersion parameter (vq = 72 km/s) than the high velocity component 
(// = 270 km/s, vq = 30 km/s). To check whether these correspond to the thick and thin 
disks populations, we take subsets of the previous selected star particles with < < 
150 km/s and > 275 km/s respectively, and calculate the disc scale height as well as the 
total velocity dispersion. We obtain z^J^^'^^ ~ 3 kpc, z^'^*" ~ 0.67 kpc, a^^l'^'' ~ 165.5 km/s, 
(jthm ^ gg^g km/s. These values are higher than those measured for the Milky Way: 
from SDSS tomographic results, ^^ly — 1 kpc, — 0.35 kpc [96]; for the velocity 

dispersion, a^^w — km/s [97], ct^*^ ~ 40 km/s [98]. Although the thick and thin 
stellar disks could be revealed by applying severe cuts on Vtf,, the inferred values of the disk 
scale height and the total velocity dispersion are not reliable, they probably overestimate 
the correct values by up to a factor 2. Hence, we stay inconclusive regarding the consistency 
of stellar populations in this simulation with Milky Way observations. 

The existence of two DM components in the solar neighborhood cannot be assessed 
undoubtfully in this simulation due to the limited local resolution. Nevertheless, as shown 
on panel b of Fig. 5, a fit with a double Gaussian gives a better description of the 
distribution than with a single Gaussian. From this fit, it appears that the mean velocity 
of the rotating component is comparable for DM and for stars. The dark disk particles also 
have a smaller velocity dispersion, also, it appears that the dark disk constitutes around 
25% of the total local density in the solar neighborhood. The value of this fraction can 
be further checked by calculating average densities from the simulation. For DM particles 
in the spherical shell 7 < ii < 9 kpc, we find (pdm) = 0.31 GeV/cm^, in agreement 
with what is found in DM only simulations and the commonly adopted value in this field, 
Pdm = 0.3 GeV/cm^. However, when we restrict DM particles to a ring of thickness 6z 
around the galactic disk, the average density is higher. We find (pdm) = 0.37 GeV/cm^ 
for 6z/2 = 1 kpc and (pdai) — 0.39 GeV/cm^ for 6z/2 < 0.1 kpc. The density increase 
near the galactic disk is indeed about 25% of the total density in the solar neighborhood. 

The extent of the density increase is shown on Fig. 7. On the left panel, DM particles 
with 7 < < 9 kpc are partitioned into slices of thickness 6z = 2 kpc parallel to the 
plane xy (galactic plane), and the average density is shown as a function of z. If the halo 
was close to spherical, there would be no bump around z = 0. On the right panel, we 
show the disk to shell density ratio for different planes. For the galactic plane, the density 
increase due to the dark disk is typically between 15 and 30% for R < 25 kpc. For other 
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planes, the ratio oscillates around 0.9, except when satellite halos are encountered, this is 
the case for the plane yz and i? ~ 13 kpc, as illustrated by the density maps of Fig. 8. 
With the dark disk component, the entire DM halo appears with an oblate shape. The 
local density increase due to the dark disk is relatively mild in this simulation, compared 
to more extreme cases [71]. Therefore, following the results of Ref. [73], we do not expect 
a strong increase of the direct detection signal in this simulation compared to a standard 
Maxwellian halo. 

Finally, we checked that the DM rotation is not limited to a thin galactic plane but is 
present in a rather thick dark disk. In particular, Fig. 9 shows the variation of the mean 
tangential velocity v(j) as a function of R for 2 = (left panel), and as a function of z 
for r = 8 kpc (right), where r = \J ^ ip- is the polar radius. A significant rotation is 
found in the entire galactic plane, up to i? = 100 kpc. For r = 8 kpc, the rotation remains 
important up to heights of ~ 15 kpc from the galactic plane. 

3. Direct detection experiments and signals 

Direct detection experiments aim to measure the energy deposited during scatterings of 
Weakly Interacting Massive Particles (WIMPs) with the detector material. Given the weak 
interaction rate and the low expected signal, radiopurity of the material and control of other 
backgrounds are essential although experimentally challenging. In order to attenuate these 
difficulties, one possibility is to look for a typical signature of the signal that is hard to 
mimic with a background noise. The DAMA/Libra experiment [41] is the only one that 
has observed a positive signal with such a signature, namely the annual modulation of the 
event rate. However, as the signal could not be confirmed by any other experiment, neither 
with a heavier or a lighter target, the compatibility of the DAM A result and the exclusion 
limits set by these other experiments becomes more and more controversial [99, 100]. 

Ways of reconciliation have been sought in the detector characteristics (quenching 
factors uncertainties and channelling effects [101, 102]), in particle physics uncertainties 
(nuclear form factors [103], nucleon effective couplings [104]) or scenarios (elastic scat- 
tering [105-109], inelastic scattering [62,110-114]), and in the Dark Matter halo velocity 
distribution [43,44,62,115-117]. 

In this section, we give all the relevant steps necessary to compute the spin independent 
scattering signal of Dark Matter on nuclei in direct detection experiments. Reference values 
for all the useful parameters and quantities will be given, and discussed in light of the 
DAMA controversy. In this paper, the accent will be put on the impact of a realistic velocity 
distribution (as given by our cosmological simulation with gas and stars) on the direct 
detection signal and its modulation. In particular, we will focus on how the fit to DAMA 
and the other exclusion limits change compared to a standard Maxwellian halo. Both the 
elastic and the inelastic cases will be covered. For the influence of other parameters, the 
reader will be referred to the existing literature. 

3.1 Event rate formalism 

We consider collisions between dark matter particles from the Galactic halo and the nuclei 
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of a given low background detector. The relevant characteristics of the halo are the local 
density of dark matter, taken to have the fiducial value poM = 0.3 GeV/cm^ at the Sun's 
location, and the local distribution of velocities i{v) with respect to the Earth. 

The differential event rate of nuclear recoils as a function of the recoil energy En is 
conveniently factored as 

dU PDM da 

dE-n = M^dE-n^^'^^^'^ ^'-'^ 

where Mj^m is the WIMP mass, da/dEn encodes all the particle and nuclear physics 
factors, and t](Eji, t) is the mean inverse velocity of incoming DM particles that can deposit 
a recoil energy Er. The time dependence of the velocity distribution is induced by the 
motion of the Earth around the Sun, which leads to a seasonal modulation of the event 
rate [118,119]. 

The total recoil rate per unit detector mass in a given energy bin [Ei, E2] is obtained 
by integrating Eq. (3.1), 

nt) = J^ dEneiEit) i^—^GiER,aiER))j , (3.2) 

where e is the efficiency of the detector and the finite energy resolution of the experiment 
is taken into account by convoluting the differential rate with a Gaussian distribution with 
spread (t{Er). For detectors made of several elements, the total rate is the average of the 
rates TZi{t) for each component i, weighted by its mass fraction /j 

n{t) = Y,fMt) (3.3) 

i 

Finally the expected number of observed events per unit time is the product of the total 
rate times the detector mass Mdet- 

The particle physics is enclosed in the term da/dEn, which is generally parameterized 

as 

da_ ^ {zU + iA-Z)u] 

dER 2pl /2 

where Mat is the nucleus mass, /U„ is the reduced WIMP/neutron mass, cr^ is the zero 
momentum WIMP-neutron effective cross-section, Z and A are respectively the number of 
protons and the atomic number of the element, fp (fn) are the WIMP effective coherent 
coupling to the proton (resp. neutron), and F'^{Er) is the nuclear form factor. In the 
sequel, we will consider the case of scalar interactions, most relevant for the elastic scat- 
tering scenario, for which fp ~ /„ and the case of weak vector interactions through a Z 
boson, most relevant for the inelastic scattering scenario, for which /p = 4 sin^ 6w — 1 and 
fn = ^- The form factor F'^{Er) characterizes the loss of coherence for non zero momentum 
transfer. We use the simple parameterization given by Helm [103,120], defined as 



2 

i2 



F'{Er) , (3.4) 



F{ER) = 3e-^'^'/'^ , (3.5) 
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with q = ^J2MnEr, s = 1.0 fm, r = {(1.23(m7v/m„)i/3 _ Q_g^2 _^ q 53^2 _ 5s2|i/2 
effective nuclear radius, Ji is a spherical Bessel function of the first kind. This form factor 
is optimal for scattering on Iodine [103]. For simplicity we use the same form factor for 
all targets (more accurate form factor are off by at most 0(20%), for some targets and at 
large recoil energies [121]). 

Finally, the velocity distribution appears in the quantity 

v{En,t)= [ d^v^-^ , (3.6) 



where v the WIMP velocity wrt the Earth and Vmin is the minimum velocity needed to 
provoke a recoil inside the detector. Two cases need to be considered. In the elastic 
scattering scenario, a dark matter particle is simply scattered off a nucleus. In the inelastic 
scattering scenario, a dark matter particle DMi is supposed to scatter into a slightly heavier 
state DM2, with a mass splitting 6 = Mdm2 ~ ^dMi ~ 100 keV. In this scenario, which 
has been first proposed in Ref. [110] and confronted to the recent data in Refs. [62,111-114], 
a much broader range of dark matter candidates may both fit DAMA and be consistent 
with the other experiments. The threshold velocity is given by 



1 ( MnEr , 

with /i the WIMP-nucleus reduced mass. This formula encompasses both the elastic {5 = 
0) and the inelastic {5 7^ 0) scattering scenarios. 

When the velocity distribution in the galactic frame fgai(i') is isotropic, rj is given by 

27r 

7]=— {F{Vesc) - F{v)) dv , (3.8) 

with = min{ fgsc 1 ^min i ^©}i ^esc is the galactic escape velocity, V(^{t) is the Earth's 
velocity in the galactic frame, which varies with the time of year t, and 'P{v) = J v fgai(w) dv. 
In deriving this formula, we have used the realistic assumption Vesc > at any time t. 
Notice that F{v) is an even function of v, because fgai only depends on the modulus \v\ of 
the velocity. For a standard Maxwellian distribution with a mean velocity vq, 

fgal(tT) = — ^^^e-->0 foTV<Vesc , (3.9) 

where N(vesc) is a normalization factor, Eq. (3.8) reads 

— i— |Erf(^)-Erf(^^)-%^e-«.".=/-n , (3.10) 

2N[Vesc)Ve [ ^VqJ \VoJ VTTVo J 



which agrees with the result of Ref. [122]. It might seem curious that Eq. (3.8) is expressed 
as a definite integral between v- and f+ while particles with any velocity above Vmin 
should contribute to rj. However, the function F{v) is itself an integral on the velocity 
distribution. Also the dependence in Vesc is now explicit. In the rest of this paper, we will 
assume an escape velocity Vesc = 600 km/s, somewhat on the upper part of the expected 
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Figure 10: Variation of the phase to corresponding to a maximal Earth velocity and event rate as 
a function of the galactic circular velocity v^- 

range, 498 km/s < Vesc < 608 km/s [117]. Modifications of the velocity distributions and 
their impact on the fits, have been discussed in various works, see for instance Refs. [43, 
44,62,113]. 

The seasonal modulation of the event rate is caused by the variation of the Earth's 
velocity throughout the year. Following Refs. [57,115], in the galactic reference frame such 
that Ix is pointing towards the galactic center, ly is along the disk rotation, and Iz is 
pointing towards the galactic north pole, we write 

^^e(*) = +VEo{t) , 

Vq = Vc + Ve,pec = (0, 220, 0) + (10, 5, 7) km/s , 
VEo{t) = VEo{eismu{t - ti) - e2Cosuj{t - ti)) , (3.11) 

where the velocity vq of the Sun with respect to the halo [123] is decomposed as the sum 
of the galactic disk circular velocity Vc and a peculiar velocity VQ^pec, the Earth mean 
orbital velocity is fEo = 29.8 km/s, with u = 2tt/T and a period T = 1 year, ti = 
79.62 days is the time of the vernal equinox, and ei = (—0.067,0.4927,-0.8676), 6*2 = 
(-0.9931,-0.1170,0.01032) (the Earth's orbital eccentricity is neglected). Therefore the 
Earth velocity follows a cosine function. As veo Vq, we have 

W0(t) ~ W0 + tiEO sin7coscLi(t - to) , (3-12) 

with 7 ~ 29.3° is the angle between the ecliptic plane perpendicular axis and Vq. The 
Earth velocity is maximal for t = tQ, with 

7T 1 Vq ■ 62 

^0 = ^1 H 1 — arctan — — ~ 151.5 days . (3.13) 

2lO iO Vq ■ ei 

With these values, the Earth velocity modulation amplitude amounts to 6.2%. As pointed 
out in Ref. [124] , the circular velocity of the Sun may be larger than the value Vc = 220 km/s 
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usually considered. Also, if the DM halo is co-rotating with the galactic disk, the relative 
velocity relevant for calculating the direct detection signal is effectively reduced. In Fig. 10, 
the circular velocity is varied between 100 km/s and 300 km/s, the corresponding phase to 
increases with Vc, from to = 144 to to = 153 days. 

For an isotropic velocity distribution, the differential rate also follows a temporal cosine 
function at first order in Veo/vq- Indeed, the only time dependent factor in Eq. (3.1) is rj, 
for which we get, in the limit Vesc > v+ 

rj{ER,t)c^rjo{Ej^)+7]i{En)—sm-fcosoo{t-to) , (3.14) 

Vq 

with 

r]o{Eji) = 7]{Er,veo = 0) , 

m{En) = — r Fiv) dv - 2^(f«) + F{v^_)) , (3.15) 

with v^ = Vmin{Efi) ^ Vq. It should be stressed that the time of year t = to does not 
necessarily correspond to a maximum of the event rate. For low energy events, Eq. (3.15) 
shows that it is a minimum instead, as iji becomes negative. 

In this paper, we also need to consider discrete velocity distributions as the output of 
a simulation gives a list of particles with their position and velocity. If we denote by Npart 
the number of particles, and by Vi the velocity of the particle i in the galactic frame, the 
distribution can be written as 

fgal(^T) = ^E'^(^-^*) ' (3.16) 

-'■^part ^ 

where 6 is the Dirac delta distribution in 3 dimensions. Therefore, for 77, we get in this 
case ^ ^ 

V = — H{wi - Vmin) , (3.17) 

Npart ^ Wi 

with Wi = \vi — vq\ and H(x) is the Heaviside step function. To compute the total rate, or 
the modulation amplitude in some energy bin, we need to integrate the differential event 
rate. For discrete distributions, we take advantage of the step-like shape of the differential 
event rate for a given flow Vi to get accurate results. 

3.2 DAMA modulation data 

The former DAMA/Nal [36] and its follow-up DAMA/LIBRA [41] are experiments run 
at the LNGS at Gran Sasso, Italy using Nal(Tl) crystals as target. They are designed to 
detect the dark matter recoil off nuclei through the model independent annual modulation 
signature which is due to the motion of the Earth around the Sun. The experimental results 
obtained by DAMA/LIBRA, with an exposure of 0.53 tonxyr collected over 4 annual cycles, 
combined with the ones of DAMA/Nal, for an exposure of 0.29 tonxyr collected over 7 
annual cycles, corresponding to a total exposure of 0.82 tonxyr, show a modulated signal 
with a confidence level of 8.2 a [41]. 
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Figure 11: DAMA modulation amplitude data (in counts per day per kg per keV) as a Junction of 
the recoil energy, expressed in electron- equivalent keV units. The 36 bins, from 2 to 20 keV, show 
a modulation signal in the low energy part, from 2 to 6 ke V. The upper part is compatible with zero 
modulation, as shown in the two bins version of the data (from 2 to 6 keV and from 6 to 14 keV). 



The DAMA modulation data, reproduced on Fig. 11, are obtained by adjusting (with 
a hkelihood function) the observed event rate in each energy bin A; to a temporal cosine 
function 

Sk = So,k + Sm,kCOSUj{t-to) , (3.18) 

where Sq^^ is the constant part of the signal, Sm^k is the modulation amplitude. The DAMA 
collaboration uses the value to = 152.5 days (2nd of June) for the phase corresponding to 
a maximal Earth velocity and a maximal event rate. It should be stressed that the true 
value of to depends in general upon the velocity of the Sun, and upon the DM velocity 
distribution. However, for an isotropic distribution, Eqs. (3.13-3.14) show that the phase 
to is completely determined once the velocity of the Sun in the galaxy is specified. For the 
low energy region where the signal modulation is present, the experimental value of to has 
been determined by the DAMA collaboration. When the phase is allowed to vary freely, 
and fitting again the event rate in the energy interval 2 < Eee < 6 keV with a cosine 
function, the best-fit value for to is 

to = 144.0 ± 7.5 days (Icj) . (3.19) 

The DAMA spectrum is given in keVee (electron-equivalent keV) . The observed energy 
released in scintillation light is related to the nucleus recoil energy through a so-called 
quenching factor q, E^^^^-^ = q ■ E^^^^^^- This expresses the fact that a recoiling nucleus 
may loose energy by collisions with other nuclei, hence in the form of heat, or through 
collisions with electrons, which create scintillation light. The reference values for Iodine and 
Sodium are respectively qj = 0.09 and qjya = 0.3. However it has been pointed that the so- 
called channelled events may play a role [101, 102]. This refers to events in which a recoiled 
nucleus moves along the axis of the Nal crystal, losing most of its energy by collisions with 
electrons, in which case the quenching factor may be larger, up to g ~ 1. Once channelling 
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is taken into account, collisions of light WIMPs with Iodine become relevant, while recoils 
on Sodium become negligible in all scenarios [42,43,102,122,125,126]. We use the fraction 
/ of channelled events given in Ref. [43], 

where the recoil energy Eji is in keV. We have verified that the other parameterizations 
found in the literature [42, 112, 122] give identical results. Finally, we use the energy 
resolution of the detector given by the collaboration [41] 

a(E) 0.45 

^ = -^ + 0.0091 (3.21) 



and a detector efficiency e = 1 [41, 122]. 

To fit the observed energy spectrum of Fig. 11 with a theoretical model, we use a 
goodness-of-fit (GOF) method with a metric, for the 12 first bins of data, with width 
0.5 keVee, from 2.0 to 8.0 keVee [41]. 

/ Q Qobs\2 

X' = Y. 2 ^ (3-22) 

i=l '^i 

where Si is the predicted value of the modulation amplitude in the bin number i, Sf^^ is 
the value reported by DAMA and cij is the corresponding experimental uncertainty. As 
data above 8.0 keVee are essentially compatible with no modulation signal, including them 
in the fit enlarges the allowed regions in the parameter space. For explicit numerical 
values of 5°'"', we refer the reader to Table III of Ref. [122]. As emphasized earlier, these 
values were obtained by the DAMA collaboration with the hypothesis Iq = 152.5. To be 
consistent with their analysis and to facilitate comparisons with previous analyses, we will 
also use this fixed value of to for the determination of the allowed regions in the parameter 
space, although the true phase corresponding to the maximal event rate might be different 
for some energy bins. For two parameters {M^im and a in the elastic scenario, and 5 and 
Mum in the inelastic scenario), the metric has 10 degrees of freedom, therefore the 90% 
(99% and 99.9%) confidence level (CL) corresponds to < 16-0 (resp. < 23.2 and 
< 29.6). 

3.3 Other experiments exclusion limits 

So far all the other direct detection experiments searching for dark matter are compatible 
with null results. In this section we briefly describe the experiments that lead to the most 
constraining limits on both the elastic and the inelastic scenarios. 

In the elastic scenario, light nuclei, like Aluminum [A = 27) and Silicon [A = 29), are 
more sensitive to light WIMPs scattering Md^j ~ multi-GeV, and provide the strongest 
upper bounds on the allowed parameter space favored by DAMA. In the inelastic scenario, 
which involves heavier candidates, experiments made of heavy nuclei, like Iodine (^4 = 127), 
Xenon (A = 129) and Tungsten (A = 184) are the most constraining ones. Germanium 
(A = 73) made detectors fall in between. 
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In computing the rate of Eq. (3.2) for each experiment, we uniformly assume that the 
smah number of events seen, if any, are signals from dark matter and we use the Poisson 
statistics to find the parameter space excluded at a given confidence level (CL). Integrat- 
ing over the exposure time, the upper bounds are obtained by requiring that the total 
number of events Ntot is compatible with the number of observed events at 99.9% of CL, 
while Ref. [113] uses 99% CL. This means that we take less severe constraints, leading 
to slightly smaller excluded regions. We have checked that the exclusion curves that we 
obtain reproduce fairly well the published 90% CL upper bounds for each experiment. In 
the following we describe the main features of the experiments considered, their particular 
characteristics, and the function adopted. 

CDMS : The Cryogenic CDMS experiment at Soudan Underground Laboratory op- 
erates Ge and Si made solid-state detectors. For a heavy WIMP, the Ge data are more 
constraining: we have considered the ensemble of the three released runs, with respec- 
tively an exposure of 19.4 kg-day [127], 34 kg-day [128] after cuts and a total exposure 
of 397.8 kg-days before cuts for the third run [37]. The sensitivity to nuclear recoils is 
in the energy range between 10-100 keV while the efficiencies and the energy resolution 
of the detectors are given in [106]. For the third run, the efficiency is parameterized 
as e{ER) = 0.25 + Qm{ER - 10)/5 for 10 keV < Er < 15 keV and e{ER) = 0.30 
for Er > 15 keV. The searches on Si are more sensitive to light DM candidates. For 
the run released in Ref. [129], two 100 g Si ZIP detectors were used, and the data set 
corresponds to 65.8 live days. The energy-dependent efficiency has been modeled in 
Ref. [122] as e{ER) = 0.80 x 0.95(0.10 + 0.30(£'r - 5)/15) for 5 keV < Er < 20 keV 
and e{ER) = 0.80 x 0.95(0.40 + 0.10(^r - 20)/80) for 20 keV < Er < 100 keV. Two can- 
didate events with energies Er ~ 55 keV and Er ~ 95 keV were observed, consistent with 
zero event once expected neutron background is taken into account. For the run released 
in Ref. [128], the total exposure is 12 kg-day after cuts, and the energy range has a higher 
minimum threshold 7 keV < Er < 100 keV. 

For zero observed event and Poisson distributed data, the function reduces to 



where Ntot is the total predicted number of events in the entire range of energy, for given 
values (pi,P2, . . . ) of parameters (Pi, P21 ■ ■ ■) of the tested model. Therefore, for a 99.9% 
CL limit and two parameters, the exclusion limit is obtained by requiring that Ntot is less 
than 6.9. 

XENONIO : XENONIO is a dual-phase Xenon chamber operating at LNGS, with a 
total exposure of 316.4 kg-days and 10 candidate events in the recoil energy range between 
4.5-26.9 keV [38,39]. For the data analysis, we follow Ref. [43], using the 7 energy bins 
provided by the collaboration and a for Poisson distributed data 



(3.23) 




(3.24) 
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with Nf is the predicted number of events in the ith bin, Bi = (0.2, 0.3, 0.2, 0.8, 1.4, 1.4, 2.7) 
is the expected background, Di = (1,0,0,0,3,2,4) is the number of detected events, and 
the log term is zero if Di = 0. For a 99.9% CL hmit and two parameters, the value of 
the function must be less than 13.8. As observed in Ref. [43], using a constant nuclear 
recoil scintillation efficiency ^e// = 0.19 gives a slightly underestimated energy threshold 
of 4.5 keV. Correcting for this bias leads to a slightly less severe exclusion limit. In this 
paper however, we will keep the energy bins as given by the XENON collaboration. Finally, 
we can notice that this experiment is also sensitive to inelastic dark matter in a similar 
way as the DAM A experiment, due to the close proximity of the target nucleus masses. 

CRESST : The CRESST-I and CRESST-II experiments at LNGS are mostly sensitive 
to light and heavy WIMPs respectively. In the first data release (CRESST-I) [130], the 
prototype detector module uses AI2O3 crystals as target, with a total exposure of 1.51 
kg-days covering the energy range between 0.6 - 20 keV. The energy resolution is given by 
(j{E) = Y^(0.519 keV)-^ + (0.0408 E)"^ and the collaboration reports eleven observed events 
after cuts. Therefore, for Poisson distributed data, the number of predicted events cannot 
exceed 23 at 99.9% CL. For the second commissioned run of 2007 (CRESST-II) [40,131], 
the target is changed to CaW04. The presence of heavy nuclei of Tungsten enhances the 
sensitivity to spin-independent inelastic scatterings. The energy range for the CRESST-II 
Zora and Verena detectors is 12-100 keV, with a total exposure of 47.9 kg-days before cuts 
and an acceptance of 0.9 on Tungsten recoil. The collaboration measured seven events, 
leading to a hmit Ntot < 16 at 99.9% CL. 

Other experiments are potentially relevant, like ZEPLIN, CoGeNT and TEXONO. 
Also, the total rate observed by DAMA provides a constraint that excludes part of the 
parameter space. However, to avoid the cluttering of the figures, we do not include these 
limits in this paper, as they are weaker than the constraints presented. 

3.4 Results for the Elastic Scenario 

In the case of spin-independent elastic scatterings, the effective coherent couplings fp and 
/„ of the WIMP to the proton and the neutron have similar values. Therefore, the only 
relevant parameters are the WIMP mass M^m and the spin independent cross section on 
nucleon, here simply denoted as a (see Eq. (3.4)). The main result for the elastic scenario 
is shown on the four panels of Fig. 12, where the regions in the plane Mdm - cr favored by 
DAMA are plotted against the exclusion limits set by CDMS and XENONIO (the excluded 
regions are on the right and above the curves). From top left, the predictions for different 
halos: a) standard Maxwellian, b) this simulation, c) generalized Maxwellian with a = 1.5, 
d) Tsallis with q = 0.773. 

In the elastic scenario, two regions are compatible with DAMA. For the region on the 
right, which is totally excluded by the other experiments, the signal is due to quenched 
scatterings events on Iodine. The region on the left, due to channeled events on Iodine, 
is not totally excluded by CDMS. The recent results of XENONIO set however a stronger 
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Figure 12: Elastic Scenario , allowed regions in the plane Mdm — ct (DM mass vs. spin independent 
cross-section on nucleon at zero momentum transfer) 

a) Standard Maxwellian halo (poM = 0.3 GeV/cm'^, Vc = 220 km/s, vq = 220 km/s, Vesc = 
600 km/s) 

b) N-body simulation with baryons (pdm = 0.37 GeV/cm^, Vc = 220 km/s) 

c) Generalized Maxwellian distribution (poM = 0.37 GeV/cm^, Vc — 190 km/s, vq — 301 km/s, 
a = 1.5) 

d) Tsallis distribution (puM — 0.37 GeV/cm'^, Vc = 190 km/s, vq = 267.2 km/s, q = 0.773) 

For c) and d), the circular velocity Vc has been reduced compared to the standard value in order 
to take into account a halo rotation. DAM A contours correspond to 90, 99 and 99.9% CL. Stars 
indicate local best-fit points. All other exclusion curves are at the 99.9% CL. 



limit, so that the elastic solution becomes very marginal. Here we have chosen to fit the 
DAMA data on 12 bins, which yields a smaller allowed region than with 36 bins. For a 
standard Maxwellian halo, with vq = 220 km/s, it appears that the 99.9% CL region of 
DAMA is totally excluded by XENONIO at 99.9% CL. 
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Figure 13: Elastic Scenario Modulation amplitude as a function of recoil energy (left) or time 
(right) for the best-fit points of Fig. 12. In the case of the time variation, events with a recoil 
energy 2 < Er < 6 keV (in electron- equivalent units) were considered. 



For the simulation halo, the direct detection prediction is made by analyzing the phase- 
space of the central and best resolved Milky- Way sized galactic DM halo. The local struc- 
tm'e is obtained by selecting the particles {Nring = 2, 662) located in a ring (7 < i? < 9 kpc 
and |2;| < 1 kpc) around i?o = 8 kpc in the galactic plane. Despite the small number of 
flows in this selection, which causes some numerical noise and artifacts in the plots, 
a few observations can be made. The DAMA regions are slightly smaller than in the 
Maxwellian case, and the best-fit points move to the left. Also, the region excluded by 
CDMS-Si slightly enlarges towards smaller cross sections. Moreover, the channeling region 
extends outside the XENON exclusion limit, so that the compatibility between DAMA and 
the other experiments is improved with the realistic halo from the simulation. While the 
best-fit point of the channeling region is still excluded, the 99.9% and 99% CL regions are 
not totally excluded anymore. 

Several factors contribute to this improvement. The departure from a strict Maxwellian 
distribution has some impact, although moderate, as shown in panels c and d of Fig. 12. 
We considered a generalized Maxwellian halo with the two values a = 1.5 and a = 1.95, 
and a Tsallis halo with q = 0.773. Their distributions are shown on Fig. 4. As in the 
simulation, a displacement of the best-fit points towards smaller masses compared to the 
standard Maxwellian case is noticeable. On the contrary, if we take a Maxwellian halo with 
parameters tuned to fit as best as possible the simulation halo {pdm = 0.37, vq = 239 km/s, 
and smaller circular velocity Vc = 190 km/s to reproduce the global average halo rotation 
velocity), the 99.9% CL DAMA region is stiU totally excluded by XENON. Other possible 
factors include the velocity dispersion anisotropy, shown on Fig. 5, and a general anisotropy 
in the velocity field, induced by the co-rotating DM component, which has a small lag 
velocity viag — 75 km/s compared to the galactic disk. A detailed modeling and discussion 
of the impact of these features will be presented in a subsequent paper [132]. 
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The panels of Fig. 13 show the differential rate modulation for 2 < Eee < 6 keV and 
the time dependence of the total signal in this energy range for the two best-fit points of 
Fig. 12. Large fluctuations are seen in the modulation amplitude as a function of energy, 
they can be interpreted as numerical noise due to the relatively small number of flows used 
to describe the local structure of the simulation halo around the Sun. On the contrary, the 
integrated rate between 2 and 6 keV is a smooth, cosine-like function of time, with a peak 
around t ~ 150 as expected. The number of flows is however too small to reliably indicate 
any deviation from to = 152.5. 

3.5 Results for the Inelastic Scenario 

The dominance of inelastic interactions over elastic ones is achieved in a natural way when 
the scattering is mediated by a massive gauge boson [113]. In the simplest scenario, the 
mediator is the weak SU (2) l boson Z, so that no new gauge boson needs to be introduced 
beyond the Standard Model. The free parameters are therefore the mass M^m of the DM 
candidate, and its mass splitting 6 with the next-to-lightest DM state. In this paper, we 
will only consider this situation. More involved models, where the cross section differs from 
the weak cross section az can be found in Ref. [113]. 

Results for the inelastic scenario are summarized on Fig. 14. The DAMA favored 
regions are shown together with the exclusion limits from CDMS and CRESST-II (the 
excluded region extends on the left of each curve). As in the elastic scenario, the four 
panels are the predictions for different halos: a) standard Maxwellian, b) this simulation, 
c) generalized Maxwellian with a = 1.95, d) Tsallis with q = 0.773. 

In the inelastic scenario, for a standard Maxwellian halo {poM = 0.3 GeV/cm^, vq = 
220 km/s, 

Vesc — 600 km/s), there are three 99% CL regions compatible with DAMA, 
which are all mainly due to quenched scatterings on Iodine. The region with 5 ~ 30 keV 
is excluded by the constraints from both CDMS and CRESST. It is actually excluded 
by DAMA itself, as it leads to a year averaged rate higher than observed. The region 
with 5 > 50 keV and Mom > 1 TeV is the largest in terms of parameter space. For a 
standard Maxwellian halo, this region is only marginally compatible with the exclusion 
limits of CDMS and CRESST. Let us notice that in a minimal framework with only the 
two DM states, and their coupling to the Z, the standard out-of-equilibrium freeze-out 
picture leads to a DM relic abundance higher than the value measured by WMAP [133], 
because the (co) annihilation rates are suppressed for Mom ^ Mz- For scalar candidates, 
the adjustable mass of the charged SU{2) partner of the DM candidate enables to obtain 
the correct relic density for a whole range of mass at the TeV scale [134]. Finally, the 
DAMA region with Mom < 100 GeV is not excluded by other experiments. For these 
candidates, due to the proximity of the Z pole, the relic abundance is lower than needed 
for WMAP unless some asymmetry in the dark sector prevents DM density from collapsing 
during the freeze-out [45] . 

The impact of the numerical noise from the finite number of flows in the simulation 
becomes more severe in the case of the inelastic scenario. Indeed, for inelastic interactions, 
only particles with a velocity high enough can produce a recoil. From Eq. (3.7), it follows 
that the minimum velocity i)* corresponding to the energy E^, = pd/M^ is given by v.^ = 
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Figure 14: Inelastic Scenario , allowed regions in the plane 5 — Mom (DM mass splitting vs. DM 
mass) 

a) Standard Maxwellian halo (pdm = 0.3 GeV/cm'^, Vc = 220 km/s, vq = 220 km/s, Vesc = 
600 km/s) 

b) N-body simulation with baryons (pdm = 0.37 GeV/cm^, Vc = 220 km/s) 

c) Generalized Maxwellian distribution (poM = 0.37 GeV/cm^, — 190 km/s, vq — 332 km/s, 
a = 1.95; 

rfj Tsallis distribution (poM = 0.37 GeV/cm'^, Uc = 190 km/s, vq = 267.2 km/s, q = 0.773j 
For c) and d), the circular velocity Vc has been reduced compared to the standard value in order 
to take into account a halo rotation. DAM A contours correspond to 90, 99 and 99.9% CL. Stars 
indicate local best-fit points. All other exclusion curves are at the 99.9% CL. 



V^25//i. As a consequence, in the inelastic case, the signal is due to the high velocity tail 
of the distribution. As panel b of Fig. 14 shows, the DAMA region with 5 ~ 30 keV is 
strongly affected by numerical noise. In this region, only a small number of flows with a 
velocity around ~ 235 km/s contribute to the time dependent modulation. The region 
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Figure 15: Inelastic Scenario Modulation amplitude as a function of recoil energy (left) or time 
(right) for the best-fit points of Fig. I4. In the case of the time variation, events with a recoil energy 
2 < Eft. < 6 keV (in electron- equivalent units) were considered. 



with Mdm < 100 GeV is also reduced, as it requires flows with a velocity higher than ~ 
700 km/s. The third region is probably the most reliable. Its deformed shape, compared 
to a Maxwellian halo, as well as its relative position to the exclusion curves are features 
that show that the velocity distribution of the simulation halo is not Maxwellian and not 
isotropic. In particular, the compatibility between DAMA and the other experiments is 
strongly improved. DAMA solutions at 90% CL are found, while only solutions at 99% 
CL were available in the case of the standard Maxwellian halo. This improvement of the 
compatibility can also be attested by the different position of the best-fit point for each 
case. 

Results for a generalized Maxwellian or a Tsallis halo are shown on panels c and d 
of Fig. 14. The deviations from Maxwellianity cause the DAMA regions to shrink. Most 
importantly, in the region 5 ~ 120 keV and Mum — 1 TeV, the strong improvement of the 
fit between DAMA and the other experiments seen in the simulation is reproduced. For 
this region, the relative position of the best-fit point and the exclusion curves of CDMS 
and CRESST-II is comparable for the Tsallis and the generalized distributions, as both 
give a good fit of the high velocity tail of the velocity module distribution found in this 
simulation (see Fig. 4). However, a Tsallis distribution favors slightly larger mass splittings 
6, and can be seen to be in better agreement with the simulation. 

As in the elastic case, the differential rate modulation for 2 < Eee < 6 keV and the 
time dependence of the total signal in this energy range for the three best-fit points of 
Fig. 14 are shown on Fig. 15. The impact of numerical noise is clearly seen in the time 
modulation, which can strongly differ from a cosine-like behavior. 
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4. Summary &; Perspectives 



In this paper, we have analyzed the phase-space structure of a galactic halo extracted from 
an advanced cosmological N-body simulation which contains stars, gas, and DM, and used 
this information to make direct detection predictions. 

Usually, such predictions are done with simplified astrophysical assumptions. The local 
phase-space structure in the solar neighborhood is taken as a smooth, isotropic Maxwellian 
distribution. However, it is known for quite some time that deviations should be expected 
in any realistic DM halo. The action of long-range gravitational forces, as well as the 
presence of a non-virialized halo component in the form of cold streams from the contin- 
uous secondary galactic infall [115, 116] already contribute to distort the simple picture 
of an isotropic Maxwellian halo. Cosmological numerical simulations done in the ACDM 
paradigm have demonstrated the hierarchical character of structure formation, leading 
to galactic structures with numerous sub-halos [135, 136]. Recent N-body simulations 
with very refined resolution, such as Aquarius [68] or Via Lactea [67] have even shown 
that non-Gaussianity and anisotropy (the velocities are rather in the radial direction) are 
present [43,54,62,137]. However, these simulations contain only DM particles, and cannot 
therefore provide a realistic description of a galaxy in a satisfactory way, despite the large 
number of particles. Dynamical interactions between the baryonic and the dark compo- 
nents of a galaxy are indeed expected to cause non negligible distortions on the phase-space 
structure of both components [86] . Recently, simulations with baryons have shown that the 
presence of a galactic disk drags merging satellites in a preferential direction [71], which 
leads to the formation of a thick stellar disc and a so-called dark disk. 

In this simulation, it appears that the DM halo contains a dark disk component that 
is co-rotating with the galactic disk. The local DM density in the solar neighborhood has a 
value pdm = 0.37 to 0.39 GeV/cm^, consistent with the recent determination of Ref. [49], 
but higher than the usually quoted value of 0.3 GeV/cm'^. Despite large fluctuations, the 
average circular velocity differs significantly from zero throughout a thick region of the 
halo. At the Sun's location, a double Gaussian provides a reasonable fit of the distribution 
of the tangential velocity u^. We infer that the rotating DM component has a mean lag 
velocity viag — 75 km/s compared to the stellar disc, and contributes to around 25% of the 
total local density. This rather small value, compared to results of other ACDM simula- 
tions [71], rather points towards a quiescent merger history for the galactic halo considered 
in our simulation. For the other velocity components and for the velocity module, we also 
observe strong deviations from Gaussian and Maxwellian distributions, which we parame- 
terized by introducing different generalizations of Gaussian and Maxwellian distributions, 
see Eqs. (2.1-2.3). For the velocity module of DM particles around = d> kpc, we find 
that a Tsallis distribution with a parameter q = 0.773 gives an excellent fit. 

The direct detection predictions are given as fits for the DAMA modulation signal and 
exclusion limits for null experiments. Both the elastic and the inelastic scenarios have been 
considered with spin-independent cross-sections. As opposed to previous works [43,62], here 
the regions in the parameter space are obtained directly from the velocity distributions 
extracted from the simulation data, without modeling. They are then compared to the 
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predictions of various analytical models. 

As the co-rotating dark disk contribution to the local density is rather small, the ab- 
solute detection rates are not significantly enhanced compared to the standard Maxwellian 
case. Nevertheless, the compatibility of DAMA vs. null experiments shows some improve- 
ment. For the elastic scenario, the rotating DM component causes the channeling region 
to enlarge outside of the XENON constraint. However, the overall goodness of fit of a 
global solution is still poor as most of the channeling region remains excluded. For the 
inelastic scenario, the situation is much brighter, as two solutions with 6 ~ 130 keV are 
permitted. The solution with a heavy candidate {Mom ^ TeV) becomes much more 
acceptable with the realistic halo provided by our simulation. Moreover, such candidate 
can also achieve naturally the relic abundance needed for WMAP [134], and is very little 
constrained by accelerator data [138]. If the DAMA modulation signal is taken seriously, 
all these hints concur to some interesting new physics beyond the Standard Model at the 
TeV scale. The solution with a mass around Mom — 80 GeV is also acceptable from both 
the direct detection data and the particle physics point of view. A correct relic abundance 
however requires some asymmetry in the dark sector [45] or some non thermal production 
mechanism. The accelerator constraints could be more stringent, but have to be analyzed 
within a particular model [139]. These questions are beyond the scope of the present paper. 
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